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Summary 

A time domain numerical scheme is developed to solve for the unsteady flow about a flat plate 
airfoil due to imposed upstream, small amplitude, transverse velocity perturbations. The govern- 
ing equation for the resulting unsteady potential is a homogeneous, constant coefficient, convective 
wave equation. Accurate solution of the problem requires the development of approximate bound- 
ary conditions which correctly model the physics of the unsteady flow in the far field. A uniformly 
valid far field boundary condition is developed, and numerical results are presented using this con- 
dition. The stability of the scheme is discussed, and the stability restriction for the scheme is 
established as a function of the Mach number. Finally, comparisons are made with the frequency 
domain calculations by Scott and Atassi, and the relative strengths and weaknesses of each ap- 
proach are assessed. 
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1. Introduction 


In recent years, due to the availability of large, high speed computers and improvements in 
approximation methods, significant progress has been made in computational fluid dynamics in 
the effort to solve unsteady aerodynamic flow problems. Much of the work that has been done in 
this area has used the so-called primitive variable approach, wherein a system of governing partial 
differential equations such as the unsteady Euler or Navier-Stokes equations are solved in time along 
with certain specified boundary conditions. There are a number of difficulties associated with this 
approach, including uncertainties about far field boundary conditions, the need for sophisticated 
algorithms to solve the equations, and the requirements for large computer memory and CPU time. 

For many aerodynamic flows of practical interest, the unsteadiness in the flow arises due to 
small, upstream vortical disturbances (gusts) which are convected downstream and interact with 
an airfoil or cascade of airfoils. The convected vortical gusts induce unsteady pressure fluctuations 
on the blade surfaces which can cause forced vibrations and radiate noise into the far field. For 
this class of “weakly rotational” flows, it is not necessary to solve the full governing equations 
of fluid motion, but rather a linearized treatment which employs the so-called “rapid distortion” 
approximation can be used [1]. For potential mean flows, the problem can be formulated in terms 
of a single, non constant-coefficient, inhomogeneous, convective wave equation [2,3]. 

Among the advantages of the linear approach are the fact that a variety of standard differencing 
schemes are available for linear problems, the algorithms are simpler, and computational times can 
be an order of magnitude less than for the primitive variable approach. In addition, it is easier to 
derive far field boundary conditions which correctly model the physics of the unsteady flow. 

Recent articles by Atassi and Scott [4], and Scott and Atassi [5] have described a general 
numerical scheme which has been developed to implement the linear approach for the purpose of 
obtaining frequency domain solutions to small amplitude, periodic vortical flows about isolated 
airfoils. In [4] they show how the motion of a nonuniform flow about a row of moving blades 
can be modelled as a three-dimensional vortical disturbance convected by the mean flow. They 
present the general theory for the linearized treatment of such flows, and present unsteady results 
for symmetric airfoils of arbitrary thickness at various Mach numbers. In [5] they show how 
the problem simplifies if the thickness of the airfoil goes to zero, and present the details of their 
frequency domain numerical scheme for the case of thin airfoils in subsonic flows. Comparisons of 
their numerical results with known analytical solutions to the classical thin airfoil gust response 
problems indicate that their scheme can accurately predict the unsteady lift on thin airfoils for a 
large range of Mach numbers and reduced frequencies. 

The alternative to the frequency domain formulation reported in [4] and [5] is to solve the 
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linearized problem in the time domain using a time marching procedure. There are some possible 
advantages to the time domain approach. For example, in the frequency domain approach, it is 
necessary to solve a large matrix equation using a direct solver. Due to the time and storage 
requirements of direct solvers, a time marching scheme could be potentially more storage and 
computationally efficient. In addition, the frequency domain approach is limited to solving the 
problem for a single frequency and then superposing solutions for more general disturbances in 
which multiple frequencies are present. For vortical disturbances that contain a large number of 
harmonics, many calculations would be necessary and the frequency domain approach would be 
inefficient. The time domain approach, however, can handle the multiple frequency case in a single 
calculation, and would in this case be more efficient. 

The purpose of the present paper is to present the details of a time domain numerical scheme 
which has been developed to solve the linearized vortical flow problems as specialized to the case 
of thin airfoils in a transverse gust. In Section 2 we derive the linearized equations and formulate 
the boundary value problem for the case of flat plate airfoils. In Section 3 the far field boundary 
condition is derived, and in Section 4 the basic numerical scheme is presented. Section 5 discusses 
the details of the calculation of the unsteady response function, and in Section 6 we discuss com- 
puted numerical results, and compare with the results of Scott and Atassi as reported in [5] . 

2. Formulation of the Problem 

2.1 Derivation of Governing Equation 

The formulation that is presented here follows that given in the recent review paper of Atassi [1], 
Consider the equations of fluid motion for an ideal gas which is inviscid and nonheat conducting. 
The governing continuity, momentum, and energy equations may be written 


Pt + V • (pv) = 0 

(1) 

v* + (t? • V)v + ivp = 0 

O 

( 2 ) 

r 

si + t? ■ V« = 0 

( 3 ) 


where p, v , p, s are the fluid density, velocity, pressure and entropy respectively. 

This is a system of highly nonlinear equations. Since we assume that the upstream vortical 
disturbances are small compared to the mean velocity, we linearize about the mean flow state and 
introduce perturbation quantities as follows. Let 

p = po + ep (4) 

v = vo + etZ (5) 
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P = Po + ep 
8 — so + es 


( 6 ) 

( 7 ) 


where e is small compared to the corresponding zeroth order quantities that govern the mean flow 
whose density, velocity, pressure, and entropy are po, vo, po, and so> respectively. 

In the case of thin airfoils at zero mean incidence, the mean velocity is a uniform parallel flow 
in the x direction, so that vq — voe x . The first order asymptotic expansion ( O(e) ) is thus 



?j- 

si «■ 

“1? 

II 

1 

< 

(8) 


,d d w 

Po< ~m +Vo r I )u = ~ Vp 

(9) 


,d d % 

(10) 


( at +vo fe )5 = °- 

where cq — \/7Po/po 

is the acoustic speed. 


Following the traditional splitting of the unsteady velocity into solenoidal and irrotational 

com- 

ponents [1], let 

3 = ui + V<j> 

(ii) 

where t?i satisfies 

V • 3 1=0 

(12) 


and 


,d d _ 

( ai + ,, "ai ) “ 1 = 0 - 

Using equations (11) and (13) in (9), there results 


(13) 


d d 


(14) 


M^ + vo^H = -p + f(t), 

where /(f) is an arbitrary function of t. At upstream infinity <f>, <f> x and p vanish, 
deduce that 

,d d . J 

+ v °di^ = 

Substituting (15) into (8) and using (11) one gets 


From this we 
(15) 




(16) 


Equation (16) is the governing equation for the unsteady potential for the problem of small distur- 
bances to a uniform flow past a flat plate airfoil. 
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The total velocity field t? is given by 

t? = VQ?x + Ui + V<£. (17) 

We observe that once the potential field is determined, the pressure and the velocity are determined 
from equations (15) and (17), respectively. 

2.2 Boundary Conditions 

Now consider the problem of flow past a flat plate airfoil at zero degrees mean incidence with 
small amplitude, unsteady velocity disturbances imposed upstream (See Figure 1). The mean flow 
is a uniform parallel flow in the x direction, and the y direction is taken to be perpendicular to the 
airfoil. Now from conditions (12) and (13) we conclude that fli must be of the form 

H 1 = Ui(x-v 0 t,y) (18) 

and that the upstream disturbances are convected by the mean flow. Here 

fli = we x + ve y , (19) 

where w and v must be chosen to satisfy the divergence free condition (12). The unsteady velocity 
u must satisfy u — * tti as x — ► — oo. Therefore <f> must satisfy 

V<f> —* 0 as x — * — oo. (20) 

At the airfoil surface, the normal component of the velocity must vanish. This leads to the 
requirement that 

v + —^ = 0 . ( 21 ) 

ay 

In the wake behind the airfoil <j> is not continuous but must satisfy a jump condition determined 
by the continuity of the unsteady pressure. Applying (15) on each side of the vortex sheet behind 
the airfoil leads to the condition 

(^+t'°^)A^ = °, (22) 

where A <f> is the jump in <f> across the wake. 

2.3 Nondimensionalization of the Problem 

For computational purposes it is convenient to nondimensionalize the problem. We therefore 
introduce the following normalization. Let 

T = m* ( 23 ) 
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* ~ 1/2 


Y = 


_V_ 

1/2 

fo 

v 0 . 


(24) 

(25) 

(26) 


In the above nondimensionalization, / is the chord length of the airfoil. Denoting T, X, Y by t, x 
and y again respectively, the problem takes the following form: 


at 2 


+ 2 

d_ 

dt 


dtdx dx 1 

c*VV 

governing equation 

(27) 


0 

in the wake 

(28) 

d . 

= 

0 

airfoil surface 

(29) 

V<t> -> 

0 

as x — ► — oo. 

(30) 


Note that the quantity c is now the inverse Mach number of the mean flow. 


3. Absorbing Boundary Conditions 

The boundary value problem prescribed in equations (27)- (30) is an open domain problem. In 
order to obtain numerical solutions to this problem it is necessary to truncate the domain onto a 
finite region as shown in Figure 2. Because we consider only flat plate airfoils with zero thickness, 
the solution <f> is an odd function with respect to y. Therefore it is only necessary to solve the 
problem in the upper half plane. The computational domain is the rectangular region defined by 


-L, L] X [0, H]. 


Because the open domain problem has been truncated onto a finite region with fixed boundaries, 
and since the gust response problem is essentially a wave propagation problem, it is necessary to 
derive far field boundary conditions which allow outgoing acoustic waves to pass through the far 
field boundary without being reflected back into the computational domain. Since the governing 
equation is a convective wave equation, and not the regular wave equation, there are difficulties in 
deriving nonreflecting far field boundary conditions. We shall present an approach due to [8] which 
is used in our computations, and we will now describe this approach in detail. 

Generalization of the ideas to follow are found in Hariharan and Hagstrom [8]. This procedure 
has similarities to that of Bayliss and Turkel [6], but is obtained through a systematic derivation 
and has general validity. The condition which is developed in [6] is valid only in the outflow, 
whereas the conditions developed in [8] are valid uniformly in all directions. 

To present the ideas, we consider the convective wave equation in its original form 


<f>tt + 2vq <f>zt + Vq4> xx = CqV V. 


(31) 
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The idea has both a geometrical and an algebraic point of view. The geometrical view is that we 
seek a transformation in which the waves at infinity are cylindrical with respect to a moving source 
convected along the x axis with a speed U ( ». We shall now describe this process. 

Suppose we consider a wave front at a time t = t*. From time t = 0 to this time a source in 
the flow is convected a distance vot* . The signals would have propagated to a distance cot*. Thus 
the first step in deriving the transformation we seek is to develop a relation between vot* and c^t* 
with respect to a fixed observer. This is depicted in Figure 3. 

In this figure, if we label the distance from the origin to the wave front by Rcot*, we obtain 

CqR 1 + v 2 — 2Rvo cos 0 = c 2 (32) 

Solving for R we obtain: 

R{9) = Moo cos 9 ± y'l-M^sin 2 * (33) 


where M , We conclude on the basis of the behavior near 0 = 0 and 9 = •k that only the 
positive sign is valid. 

Thus the desired transformation in which the cylindrical properties are obtained is 


c 0 R(9)’ 


A 

9 = 9. 


(34) 


The interpretation is that at r = co R(9)t*, f = t* =constant, the wave front is cylindrical. The 
transformation appropriate to the seeded equation 


4>tt + 2 <f>xt + 4>xx — c 2 V 2 ^ 


(35) 


is 


where 


S{9) = 


r = rS(9 ) 
1 


(36) 


(37) 


cos 9 + %/c 2 — sin 2 9 

Now in the new coordinate system (f, 0), which we call the wave front coordinate system, one 
can represent the solution in the far field following Friedlander [7] in an asymptotic form. That is 


+ ■ 


fit - f) 
fl/2 


(38) 


From this, we obtain a radiation boundary condition that is similar to the one that is reported in 
[1], which is 

<f>t + tr+^=0 to O(rt). (39) 
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Thus the desired boundary condition (which we call the H-H condition) is a translation of (39) in 
the original coordinate system. This results in the boundary condition 

l 4> 

*• + S(S)*' + MW = °- 
This condition is implemented in the results reported in the present paper. In the numerical 
calculations this is implemented as 

<f>t + ;^yl cos W* + sin (^v ] + = 0. (41) 

We remark here that, in contrast to the geometric derivation given above, the transformation 
(37) can also be obtained algebraicly by solving the Eikonol equation corresponding to (35) which 
is 

1-2 <r x +al = c 2 (<rf + cr 2 ). (42) 



4. Numerical Scheme and Stability 


4.1 Basic Scheme and the Stability Constraint 

The indices used in our numerical calculations are » and j, where * varies from 0 to M along 
the x axis and j varies from 0 to N along the y axis. The index n is used to denote values at the 
time nAi, where At is the time step. The grid spacing in the x and y directions is denoted by Ax 
and Ay, respectively. 

Unlike the regular wave equation, the convective wave equation yields difficulties in obtaining 
accurate schemes due to the presence of the mixed derivative term <f> x t . A scheme that is second 
order accurate in space and in time was obtained using the following approximations: 


ttvvKi 


*8 l - Wi + * 




At 2 

A"+l _ AJ>n+ 1 


WT - W 

4AtAx 

Ax 2 

6 " - 2 d>?. + 6?. , 

Ay 2 


(43) 

(44) 

(45) 

(46) 


Note that all the second order derivatives are evaluated using the standard central difference formula 
except the mixed derivative term. The central difference formula cannot be used for the mixed 
derivative term, as it results in an unstable scheme. The truncation error corresponding to this 
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term can be shown to be equal to 


Af 2 At 2 

ERRORtr une ation = -—^Vtttt — 

Ay 2 


2Ax 2 


Ax 2 


<f>txxx + ~75“( c - l)4>xxxx + ~^ 2 ~ c2< f > \ 


3 12 

= 0(At 2 ,Ax 2 ). 


rvvw 


( 47 ) 


Substituting equations (43) - (46) into (27), we obtain the following numerical scheme: 

C'-^+C 1 

+0.5B(^t> - 4^+u + - 3C‘ + - *tlj) (48) 

+(1 - c J )fi 2 (^" +1J . - + ti-ij) - - 2^”j- + j) = 0, 


where 

_ At _ At 
A* Ay' 

To discuss the stability of the scheme we let the error term be the form 

Efj = A^e iaAxl e ipAvj . (49) 

Substituting into (48), we get the error equation 

£ 2 (1 + 1.5# - 2Re' iaAx + 0.5Re~ i2aAx ) 

+ e [-2 + 4(c 2 -l)# 2 sin 2 (^) + 4c 2 R 2 sin 2 (^)] 

+ (1- 1.5# + 2#e-’“ Ax - 0.5#e-‘' 2 “ Al ) = 0. (50) 


We solve this equation numerically for all possible aAx and /?Ay for fixed R and c, and find that 
the roots of the equation are all within the unit circle if the following stability criteria are satisfied: 
for Moo = 0.8, # < 0.15, 
for Moo = 0.5, # < 0.12. 

Note that there are fictitious points, namely those points at j = —1, j = N + 1, * = -2,-1, 
and * = M + 1. Thus, the scheme given by equation (48) is valid only in the interior domain, and 
special treatment is required on the grid line corresponding to i = 1 near the upstream boundary 

(«• = o). 


9 


4.2 Numerical Differencing Scheme 

The finite difference equations used to model the governing equation and various boundary 
conditions are presented below. 


4.2a Upstream Boundary 

The calculation is begun at the left boundary, where the boundary condition is 

* + s^y 1 C °’ H ‘ + ainH ' 1 1 + 2^m = °- 

The differencing used at the left boundary is 

(*)o J = ^(C 1 - tfW 

(«„• = -Wj + 4 #j - th ) 

(<t>v)oj = toj+i - 4>oj- 1)» 

where 


(r)oj = yjx l + yj 


(cos{0)) 0 j = ^2- 

r 0,3 

(sin(<?)) 0 j = ~ 

Tr\ J 
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W) "Sbj 


f 0j 


+ \/ c2 - y i/ r oj- 


The finite difference equation for the boundary « = 0 is then given by 


<' = TT 1 !!- < 

1 + ^7 

~n~q~ t *o(-3^oj + 4 <£lj ~ <t>2,j) 

Zr OJ b OJ 

+yy(^oj+i - rffc-i) + ^o,7 1 ] } i = l to jv - 1. 


(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 

(58) 

(59) 

(60) 


4.2b Interior Points 

For interior points the governing equation is 


<f>tt + 2 <f> x t + <t> zx = c 2 V 2 ^. 


(61) 
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We use an 0(At 2 ,Ax) scheme for those points with * = 1, for j — 1 to j = N. The difference 
approximations are 


(Ml = ^5 WJ 1 - ^*1 + «;') (62) 

(Mid = ^1 - **1 + «j) (63) 

(Ml = + 4>l-i) (64) 

(«U = 2AIaJ WJ‘ - «’ - C + O (66) 

Therefore the difference equation for i = 1, j = 1 to j = N — 1 is 

tfj 1 = rh < ~ + O 

+(c 3 - l)fi 3 (^j J- - 20" j- + 00 J ) 

+ tfj-l) } (66) 


For interior points where * = 2 to » = M — 1, j = 1 to j = N — 1, the 0(At 2 , Ax 2 ) finite 
difference equation corresponding to equation (59) is 


dpi 1 = z [2 dp.- dpZ 1 

9, ‘ } 1 + 1.5 R [ 9, ’ } 9 '-’ 


h r 


- 2 (-C« + tf-i- - 3 C + 4 Ci,i - *T&) 

+(c 2 -l)f? 2 (^i >J -2< y + ei.y) 


(67) 


4.2c Airfoil Surface 

At the surface of the airfoil the boundary condition is 

d<f> n 
v + -£- = 0. 

dy 


(68) 


The corresponding difference equation is 

<t>to l = l ( - W.t 1 + KoAy ), i = M x to M 2 ,j = 0, (69) 

where Mi is the » index of the grid point at the airfoil leading edge and where M 2 is the i index of 
the grid point at the trailing edge. 
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( 70 ) 


4.2d Top Boundary 

On the top boundary the far field boundary condition is 

h + + **-'*• ) + = 0 
where the difference approximations are 


and where 


(MV 1 = ^(O* - Kn + O- 2 - O- 2 ) 

(MVi = «7r:(^f+i,jv-i - 


2Ax 

1 


(^)a-i = ^(O 1 - O- 2 + o 1 - 0 - 2 ) 


n _l yfl-1 


wa-i^o'+o 1 -*) 


( r )«. at- 1 = V*r+y^I7- 


(71) 

(72) 

(73) 

(74) 

(75) 


The difference equation is then 
tf*+i - 


<4?+* = ~ 


1 


f J.n-1 1 o.n+1 J.n—1 

, fii/w-i - At ~ ( ^f,AT + ri,N—2 ~ 9 x,N- 2 

r i,JV — l &i t N — 1 4»'i,Ar-l'Si,JV-l 


f? 


r « > Af-l-S'«,N-l 


[2*t ($+ 1,AT— 1 “ ^f-l.JV-l) 


Ax 


+iW-l(“0-2 + O 1 “ O- 2 ) + “0-2 ] } 

for t = 1 to M — l,j = 1\T. 


(76) 


4.2e Downstream Boundary 

At the downstream boundary the far field condition is 


+ Tm 1 + > + M(ii = 0 


where the differencing used is 


(77) 


- «rj + KUj - «-*«) 

(«Jr-w = 4^ Wi - «?-*« + «ri - «-*«) 




{4>v)m- l,j 
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(78) 

(79) 

(80) 



and where 



2 + ^Af-l,j) 

(81) 

(r)M-ij = 


(82) 


The finite difference approximation 


id biigu 


( ~ ^M-2,/ + ^M-2,; 


±r i+l __ £ 

j , -Rzm-i j At 


Arc 

+ 2 J/;($vf-i,;+i ~ 1 } 

for * = M, j = 1 to N — 1. 


(83) 


4,2f Wake Boundary and Corner Points 

For the wake boundary condition, we make use of the fact that <f> is an odd function of y, so that 
A <f> = <f>+ - <f>~ = 2^ + , where + and — denote quantities above and below the wake, respectively. 
The wake boundary condition (28) then becomes 


d . d . ^ 

at 9 ^ ax 9 


(84) 


and the difference equation is 


C‘ = +% l - m+,.o - «... j 


for » = M 2 to M — l,j = 0. 


(85) 


For the last grid point in the wake (i = M, j = 0), we use 


j.n+1 _ 


1 + ii 


[ + ^M-1,0 - ^Af-1,0 + ^(^M-1,0 + 1,0 “ ^M,o) ] 


n+1 


( 86 ) 


Lastly, for the corner points 4>\f t Ny <f>o,N we just take <f> to be the average value of the two nearby 
points, i.e., 


^ | V = 0 . 5(^+ 1 _ 1 + ^+ 1 ) ( 87 ) 

= + ^M-l , n )' (88) 
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5. Calculation of the Unsteady Response Function 


For purposes of comparing numerical results with known solutions to the classical thin airfoil, 
gust response problems, we introduce (following Sears [9]) the unsteady response function 

R{ki, Moo) = rplvQe {iut) > ( 89 ) 

where L is the unsteady lift, and k\ = where u is the angular frequency of the incident 
disturbance, is the reduced frequency. The lift is calculated by integrating the unsteady pressure 
over the surface of the airfoil. 

Since we solve the problem in the time domain, and not the frequency domain, we obtain the 
response function through the use of the Fourier transform. The upstream velocity disturbance is 
taken to be a wide band Gaussian pulse, and we obtain the airfoil response as a function of the 
reduced frequency as explained below. 

Suppose that the value of the disturbance v(x,y,t) at the center of the airfoil is f(t). Then the 
time spectrum of this disturbance will be 



(90) 

The response of the airfoil to the disturbance will be 


=/<:>, 

(91) 

7C PCVq 

where the lift L(t) is 


il 

** 

(92) 

By taking the Fourier transformation of g(t), we get the representation of g(t) 
domain G(w), and the frequency response function of the airfoil is 

in the frequency 

R(( .\ _ G ( w ) 
F(w)- 

(93) 

For example, let us consider the case where 


f(t) = y/2ae(~ at *h 

(94) 

The spectrum of /(t) is then 


II 

(95) 


and the pulse function (disturbance) for the problem under consideration at time t and at a point 
x is f(t — x). 
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The numerical integration involved in evaluating the transforms F(u>) and G(u>) depends on the 
value of the amplitude a. For example, for a value of a = 4, Figure 4 shows the behavior of f(t). 
Since the decay rate of this function is fast, one may choose a time interval of the form [-To, To]. 
A typical value of To for our grids was 6.975. The corresponding transformation is given in Figure 

5. For numerical purposes, the index for the time begins at t = —To so that n = 0 corresponds to 
t = —To. At the time level n=4000, the pulse has passed the airfoil and the response function of 
g(t) tends to zero. Therefore this is a good place to truncate the response function g(t) numeri- 
cally. We can then calculate its Fourier transformation G(u>) very easily using Riemann sums. The 
response function g(t) and its transforms are given in Figures 6 and 7, respectively. Finally, using 
(93), we get the frequency response function R(w). The various frequency response functions given 
in Figures 8-13 are obtained by plotting Re[R(a;)] vs Jm[f?(w)]. 

6. Results and Discussion 

6.1 Comparison of Numerical Results with Known Solutions 

In order to assess the accuracy of the numerical scheme that has been developed, we compare 
our computed results with known solutions to the classical thin airfoil gust response problems. 
Because the stability constraint of our scheme depends on the Mach number, we are limited to 
calculating cases for which the Mach number is not small. For comparison purposes we present 
results for Mach numbers of .5 and .8. 

Figure 8 shows a comparison between the classical results obtained from a Possio [10] solver and 
numerical results obtained from the present scheme for a .5 Mach number. The grid dimensions 
for the results shown in Figure 8 were 15 units in the streamwise direction and 10 units in the 
normal direction. The spacing in the x and y directions was uniform with Ai = Ay. There were 
151 points in the x direction and 101 points in the y direction. 

In Figure 9 we make the same comparison as in Figure 8, except in this case the grid dimensions 
have been extended to 22.5 by 15. The number of grid points in each coordinate direction is the 
same as in Figure 8, so that the spacing is slightly more coarse on this grid. The results shown in 
Figure 10 were run on a still larger grid, with dimensions of 30 by 22.5, and 201 points in the x 
direction and 151 points in the y direction. 

By comparing the results shown in Figures 8-10, it is seen that the agreement for the high 
frequencies is very good in all three cases. However, the agreement for reduced frequencies below 
about .3 is only fair. Note that the amount of error in the low frequency cases diminishes as the 
grid dimensions become larger. This suggests that a larger grid is required to adequately model 
the long wavelength disturbances. 
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In Figures 11-13 we present results for the case of an .8 Mach number. The grids used for these 
figures correspond to the grids used in Figures 8-10, respectively. The agreement for the .8 Mach 
number case is analogous to the agreement for the .5 Mach number case, and again we see that the 
accuracy for lower frequencies is considerably improved on the larger grids. 

6.2 Comparison of the Present Scheme with the Frequency Domain Approach of Scott 
and Atassi 

As mentioned in the introduction, Scott and Atassi [5] have reported a general frequency domain 
approach that was developed to solve periodic vortical flows around isolated airfoils. Their scheme 
has wide applicability to many different flow configurations, and was recently extended to handle 
the difficult problem of vortical flows around lifting airfoils in subsonic flows [11]. Our purpose 
in the present section is to compare the present time domain calculation procedure with their 
frequency domain approach as specialized to the case of flat plate airfoils. 

There are two main disadvantages to the frequency domain procedure reported in [5]. Perhaps 
the biggest disadvantage is that it requires the solution of a large matrix equation by a direct 
solver. Because the sparse solver described in [5] has been optimized for storage efficiency, it is not 
particularly computationally efficient. Large solution times can be required, on the order of 200 
seconds per frequency on a Cray X-MP, as the reduced frequency becomes larger and the number of 
grid points increases. The other main disadvantage of the frequency domain scheme described in [5] 
is that there is a loss of accuracy as the Mach number and reduced frequency become large. Due to 
this loss of accuracy and the longer solution times for the higher frequencies, the frequency domain 
numerical scheme described by Scott and Atassi is not suitable for the very high frequencies, i.e., 
for reduced frequencies above about 5.0. 

There are, however, a number of advantages to their frequency domain approach. It is very 
accurate for a large range of frequencies, including the very lowest frequencies. Their scheme is 
valid for both incompressible and compressible flows, and can handle three-dimensional gusts. In 
addition, the unsteady grid and far field boundary condition described in [5] are ideally suited for 
acoustic wave propagation, so that the acoustics can be accurately and readily determined from 
the solution. The scheme is also general and can be applied to airfoils with arbitrary geometry. 

Among the advantages of the time domain scheme presented in this paper are that the scheme is 
very accurate at the higher frequencies, and it can handle nonperiodic vortical disturbances as easily 
as periodic ones. In addition, the present numerical approach is both storage and computationally 
efficient. For the special case of a flat plate airfoil in a transverse gust, the scheme is very efficient, 
as the airfoil response to a continuous spectrum of frequencies can be calculated simultaneously. 
Solution times for the entire response function range from 86 to 343 seconds on a Cray X-MP, 
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versus 600 to 1500 seconds for the same calculation using the frequency domain approach. 

The improved computational efficiency of the present time domain approach is due entirely to 
the fact that solutions for a complete range of frequencies can be obtained simultaneously in a 
single calculation by Fourier decomposing the results of the airfoil response to a pulse disturbance. 
This approach carries with it the disadvantage that all of the harmonics have to be captured by 
a single grid. Since the long wavelength disturbances cannot be accurately resolved on a grid 
which is suitable for the high frequency cases, there is an inherent loss of accuracy for the low 
frequency disturbances. If, as in the frequency domain approach in [5], a separate calculation were 
performed on a separate grid for each reduced frequency, then the time domain scheme is actually 
less computationally efficient than the frequency domain scheme. 

In addition to the loss of accuracy for low reduced frequencies, the other main disadvantage of 
the present time domain approach is that it is not valid for low Mach number flows. 

6.3 Conclusion 

In conclusion, the time domain numerical scheme which has been developed in the present paper 
is both storage and computationally efficient, and can accurately calculate the high frequency 
response of a flat plate airfoil in compressible flow to a transverse gust. Improvements are still 
needed in the accuracy of the calculation at the lower frequencies, and in the selection of a more 
suitable computational grid. 

In as much as nonlinear equations such as the Euler equations are usually solved using time 
marching procedures, it is hoped that the present time domain solution procedure for the linearized, 
unsteady Euler equations will illuminate some of the issues involved in time domain calculations 
for unsteady flows. 
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FIGURE 7. - FOURIER TRANSFORM OF g(t). 
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FIGURE 9. - COMPARISON BETWEEN THE NUMERICALLY COMPUTED UN- 
STEADY RESPONSE FUNCTION AND ANALYTICAL RESULTS FROM A 
POSSIO SOLVER. 
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FIGURE 8. - COMPARISON BETWEEN THE NUMERICALLY COMPUTED UN 
STEADY RESPONSE FUNCTION AND ANALYTICAL RESULTS FROM A 
POSSIO SOLVER. 
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FIGURE 10. - COMPARISON BETWEEN THE NUMERICALLY COMPUTED UN- 
STEADY RESPONSE FUNCTION AND ANALYTICAL RESULTS FROM A 
POSSIO SOLVER. 
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FIGURE 11. - COMPARISON BETWEEN THE NUMERICALLY COMPUTED UN- FIGURE 12. - COMPARISON BETWEEN THE NUMERICALLY COMPUTED UN- 
STEADY RESPONSE FUNCTION AND ANALYTICAL RESULTS FROM A STEADY RESPONSE FUNCTION AND ANALYTICAL RESULTS FROM A 

POSSIO SOLVER. POSSIO SOLVER. 
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FIGURE 13. - COMPARISON BETWEEN THE NUMERICALLY COMPUTED UN- 
STEADY RESPONSE FUNCTION AND ANALYTICAL RESULTS FROM A 
POSSIO SOLVER. 
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